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Abstract 

Torreya grandis Fort, ex Lindl, a conifer species widely distributed in Southeastern China, is of high economic value 
by producing edible, nutrient seeds. However, knowledge about the genome structure and organization of this 
species is poorly understood, thereby limiting the effective use of its gene resources. Here, we report on a first 
genetic linkage map for Torreya grandis using 96 progeny randomly chosen from a half-sib family of a 
commercially cultivated variety of this species, Torreya grandis Fort, ex Lindl cv. Merrillii. The map contains 262 
molecular markers, i.e., 75 random amplified polymorphic DNAs (RAPD), 119 inter-simple sequence repeats (ISSR) 
and 62 amplified fragments length polymorphisms (AFLP), and spans a total of 7,139.9 cM, separated by 10 linkage 
groups. The linkage map was used to map quantitative trait loci (QTLs) associated with juvenile growth traits by 
functional mapping. We identified four basal diameter-related QTLs on linkage groups 1, 5 and 9; four height- 
related QTLs on linkage groups 1, 2, 5 and 8. It was observed that the genetic effects of QTLs on growth traits vary 
with age, suggesting the dynamic behavior of growth QTLs. Part of the QTLs was found to display a pleiotropic 
effect on basal diameter growth and height growth. 



Introduction 

Plant breeding has proven to be difficult because most 
economically important traits are quantitatively inher- 
ited, controlled by an interacting network of genes, each 
with a small effect, and environmental factors. This is 
especially true for forest trees, a group of species charac- 
terized by long life spans, high heterozygosity, and 
environmental heterogeneity [1,2]. It is ubiquitous that 
variation in phenotypic traits of forest trees fluctuates 
from year to year owing to varying contributions from 
genetic and environmental factors. Molecular technolo- 
gies extensively developed in the past several decades 
have provided a powerful tool to dissert phenotypic var- 
iation into its underlying genetic factors, known as 
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quantitative trait loci (QTLs) [3-5]. Genetic mapping of 
QTLs from linkage maps has become a routine avenue 
to identify and study the genetic architecture of pheno- 
typic traits [6]. 

There has been a considerable body of literature on the 
construction of genetic maps in forest trees, such as 
Eucalyptus [7], Salix [8], Betula [9], Populus [10-12], and 
Pinus [13,14]. Because it is difficult or impossible to 
obtain segregating populations from inbred lines, genetic 
mapping in forest trees should use strategies that are dif- 
ferent from those used in annual crops or model systems. 
These strategies include (1) controlled crosses made 
from two heterozygous parents, in which different types 
of markers are segregating [15-17], and (2) open-polli- 
nated progeny from a heterozygous tree, consisted of 
seeds that share the same mother but have different 
fathers [14,18]. In gymnosperms including conifers, 
cycads, ginkgo, and Gnetales, a seed consists of a diploid 
seed coat from the original parent sporophyte, a haploid 
megagametophyte developed from a megaspore produced 
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by the original sporophyte, and a diploid embryo devel- 
oped from the fertilization of the original sporophyte and 
a male sporophyte [19]. Given its haploid feature, less 
expensive dominant markers genotyped on the megaga- 
metophyte are equally informative to more expensive 
codominant markers and, therefore, the megagameto- 
phyte has been widely used as a tissue to generate mole- 
cular markers for genetic mapping. 

By taking advantage of megagametophyte-based geno- 
typing, we constructed a low-density genetic linkage map 
for an important but understudied conifer species, Torreya 
grandis Fort, ex Lindl, using dominant markers. The 
genus Torreya comprises 8 dioecious species, 5 of which 
originate in China [20]. T. grandis cv. Merrillii is the only 
commercially important ever-green variety for nut produc- 
tion. Despite its widespread cultivation in Southeastern 
China for thousands of years, however, the genetic back- 
ground of this variety has been poorly understood. In 
T. grandis and T. grandis cv. Merrillii, only several proto- 
cols have been established for molecular marker assays 
[21-24]. A limited set of markers has been used to esti- 
mate genetic relatedness among several varieties [25] and 
identify sexes in young seedlings [26,27]. 

The genetic map reported in this study is a first one for 
Torreya grandis. This map was constructed from an open- 
pollinated progeny, composed of half-sibs, produced from 
Torreya grandis Fort, ex Lindl cv. Merrillii. While markers 
were genotyped from haploid megagametophyte, phenoty- 
pic traits were measured from the seedlings of these half- 
sibs generated from the diploid embryos. Wu [18] derived 
a statistical model for mapping QTLs for diploid issues 
using megagametophyte-based markers. We used this 
model to map QTLs that control stem height and base dia- 
meter growth in juvenile seedlings. The dynamic pattern of 
QTL effects detected was studied by functional mapping, a 
statistical model for QTL mapping of developmental tra- 
jectories [28-32]. The identification of growth QTLs may 
provide scientific guidance for marker-assisted selection of 
economically important traits in Torreya grandis. 

Results 

Map construction 

The 16 primers screened from five random samples gen- 
erated a total of 109 markers, 91 of which (82%) were 
polymorphic. The number of polymorphic RAPD mar- 
kers produced per primer ranged from 2 to 11, with an 
average of 5.7. ISSR analysis using 24 primers created 
194 polymorphic loci in total on the genomic DNAs, 
with an average of 8 loci per primer. Fifteen polymorphic 
loci were amplified using primer ISSR-46, but only two 
using ISSR-52. For AFLP analysis, 14 pairs of selective 
primers generated 155 polymorphic loci, with E6/M4 
amplifying 24 polymorphic loci while E2/M7 and E4/M5 
amplifying three only, and an average of polymorphic loci 



produced by a pair of primers is 11. In total, three differ- 
ent types of markers generated 440 polymorphic loci. 

Of all the markers, 262 (75 RAPDs, 119 ISSRs and 68 
AFLPs) (59%) were identified to follow Mendelian seg- 
regation ratio according to chi-square tests at a level of 
> = 0.05. These markers were mapped to 10 linkage 
groups at an LOD > 6 with Mapmaker 3.0. The total 
map distance was 7,138 cM (Kosambi units) with an 
average distance between adjacent markers of 28.2cM 
(Figure 1). Linkage groups each comprise 7 - 71 loci 
(138 - 1898.5cM). Table 1 shows the distribution of 
markers in different linkage groups. Our long map con- 
structed in Torreya grandis cv. Merrillii confirms the 
detection of large genomes of coniferous species. 

Detection of growth QTL 

Growth curves of 510 young seedlings from the open- 
pollinated population over the first three years are illu- 
strated in Figure 2, where substantial variation in curve 
shape exists, suggesting that specific QTLs may be involved 
in developmental trajectories. Functional mapping was 
used to scan the existence and distribution of QTLs for 
growth traits by calculating genome-wide likelihood ratios 
(LR) of the full model with a QTL and the reduced model 
without QTL (Figure 3). By comparing the LR peaks with 
chromosome-wide critical thresholds determined from per- 
mutation tests, we identified several significant QTLs for 
stem diameter and height growth. As a first QTL study in 
T. grandis, we used a less conservative criterion at the chro- 
mosome level to claim the significance of QTLs. Four 
significant QTLs for stem diameter growth were detected 
in marker intervals Sll-2 - E8M8-2 on linkage group 1, 
1-77-3 - 1-36-6 in linkage group 5, 1-17-8 - 1-44-9 and 
E6M5-11 - E8M8-lin linkage group 9. Four significant 
QTLs for stem height are located between markers 1-56- 
10 - 1-46-10 on linkage group 1, 1-46-13 - E6M4-5 on link- 
age group 2, E6M7-2 - 1-34-6 on linkage group 5 and 
E6M8-6 - E7M3-12 on linkage group 8. Linkage group 1 
and 5 carry the QTLs for both diameter and height growth, 
suggesting that these segments may contribute to the 
genetic correlation between these two growth traits. 

Using the estimated growth parameters, we drew 
growth curves of each QTL identified for height and dia- 
meter traits (Figure 4A). In general, all genotypes 
increase in both height and basal diameter with age, 
although the extent of increase varies depending on the 
QTL detected. For the height QTLs detected in linkage 
groups 1, 2, and 8, favorable alleles for increasing stem 
height are derived from the egg of T. grandis cv. Merrillii, 
whereas the favorable allele at the QTL detected in link- 
age group 5 is from a pollen which may be generated by 
either T. grandis cv. Merrillii or another tree in the pool. 
For the diameter QTLs on linkage groups 1 and 5, the 
egg of T. grandis cv. Merrillii contributes favorable alleles 
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Figure 1 Genetic linkage map constructed from 75 RAPD, 
maternal parent in Torreya grandis cv. Merrillii. 



119 ISSR and 68 AFLP markers for 96 half-sib seedlings derived from a 



for increasing radial growth, whereas two QTLs on link- 
age group 9 have their favorable alleles from the pollen. 
As seen in Figure 4B, genetic effects on both diameter 
and height growth increase with age during the early 
stage of growth. 

Discussion 

Despite its economic, ornamental and ecological impor- 
tance, the genetic study of Torreya grandis has received 



little attention. For those underrepresented tree species 
like T. grandis, for which little knowledge about their 
genome structure and sequence is available, simple and 
cheap dominant markers have still been widely used for 
constructing genetic linkage maps. In this article, we 
present a first linkage map for T. grandis constructed 
from three types of dominant markers, RAPD, ISSR and 
AFLP. A total of 262 markers are clustered into 10 link- 
age groups, totaling 7,138 cM in length. Since T. grandis 



Table 1 Distribution of molecular markers on 10 linkage groups for a half-sib family of Torreya grandis cv. Merrillii. 



Linkage Group No. Markers RAPD ISSR AFLP Max. Distance Min. Distance Length Ave. Distance 
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Figure 2 Plots of stem diameter (left) and height growth (right) in the first growing seasons of young seedlings from a half-sib family 
of Torreya grandis cv. Merrillii. Part of this material, all with repeated measures at three time points, was used for functional mapping 
of growth traits. 



cv. Merrillii has a chromosome number of 2n = 2x = 22 
[33], our map only partially covers the genome of 
T. grandis. Also, given its large marker intervals, our 
map is quite sparse which should be filled by more mar- 
kers, especially more informative codominant markers. 
Nevertheless, a long genetic map identified in our study 
appears to support the big genome of a gymnosperm 
[34,35]. This map provides an additional fuel to study 
the structure, organization and evolution of conifer 
genomes. 

We used the linkage map constructed to identify sev- 
eral QTLs for height and radial growth by functional 
mapping. Although the results presented in this article 
are preliminary, our study is characterized by two fea- 
tures which are particularly important in genetic map- 
ping of forest trees. First, our mapping design is based 
on half-sibs from an open-pollinated individual tree in a 
natural stand. This design can overcome the limitation 
of genetic mapping for forest trees. Second, while many 
QTL mapping studies result from the association 
between markers and phenotypes observed at static 
points, our method makes use of growth data longitud- 
inally measured at multiple time points through func- 
tional mapping [28-32]. Since it was first proposed in 
2002, functional mapping has been used to detect 
growth QTLs in poplar [36] rice [37] and mouse [38]. 
Its tremendous advantage lies in its integration of biolo- 
gically meaningful growth equations [39] into QTL 
mapping, leading to interpretable results of QTL discov- 
eries [40]. All the QTLs identified in this study change 



their genetic effects over time, with some decreasing 
while the other increasing. 

Powerful genetic and molecular tools have made it pos- 
sible not only to detect the action of genes at the juvenile 
stage but also to detect their dynamic actions and interac- 
tion of individual loci or QTLs, which is of significance to 
seedling production in both agriculture and forestry. It has 
been well known that phenotypic formation results from 
coordinated action of genes and environmental factors. 
Based on a better understanding of gene actions, we can 
design an optimal strategy for seedling selection and pro- 
duction in the first few years in the nursery. Such a strat- 
egy is especially important for tree species of highly 
economic values, like T. grandis. 

The design used is a half-sib family with members from 
a heterozygous tree [14,18]. A more powerful design is to 
sample multiple heterozygous trees randomly from a nat- 
ural population and further collect open-pollinated seeds 
from each sampled tree [41]. This design constructs a 
two-stage hierarchical platform composed of the mater- 
nal genotypes and offspring genotypes. Wu et al. [41] 
implemented a sophisticated procedure of the EM algo- 
rithm to provide the simultaneous estimates of linkage 
disequilibria in natural populations and the QTL-marker 
linkage based on offspring genotypes. Such a design inte- 
grating population genetics and quantitative genetics is 
especially suitable for genetic mapping in T. grandis. 
More recently, Sun et al. [42] embedded genetic imprint- 
ing into this design, enabling geneticists to chart a com- 
prehensive picture of the genetic architecture of growth 
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Figure 3 The profile of log-likelihood ratios (LR) for stem diameter (A) and height growth trajectories (B) throughout the genome 
composed of 10 linkage groups in Torreya grandis cv. Merrilli. The genomic position corresponding to the peak of the curve is the 
maximum-likelihood estimate of the QTL location. The 95th percentile (indicated at horizontal lines) of the distribution of the maximum LR 
values obtained from 200 permutation tests is used as an empirical critical value to declare chromosome-wide existence of a QTL at a = 0.05. 
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Figure 4 (A) Age-dependent expression profiles of two different genotypes at the QTL detected in different linkage groups for stem 
diameter and height QTLs in Torreya grandis cv. Merrilli. Solid and broken curves represent QTL genotypes QQ and Qq, respectively. (B) 
Age-dependent genetic effects of the QTLs detected on stem diameter and height growth. The genetic effects are estimated as the differences 
of genotype 00 minus Qq. 



and morphological traits in forest trees, including 
T. grandis. 

Materials and methods 

Plant material 

The mapping population consisted of 510 half-sib progeny 
of T. grandis cv. Merrillii from a relatively isolated, open 
pollinated tree located in Changming Village, Sankou 
Town, Lin'an City, Zhejiang Province of China. The seeds 
harvested were piled for 7 days to soften and remove 



outermost arils and then stratified with sand in an open 
area to accelerate germination with two layers of plastic 
film covering on the pile top of the stratified seeds. When 
the radicle appeared, the seeds were potted in medium 
(peat: perlite: vermiculitm = 1: 1: 1), followed by partial 
megagametophyte extraction from germinating seedlings 
when seedlings were out of top medium. The white mega- 
gametophytes extracted were stored at -40°C for DNA iso- 
lation later. If megagametophytes are insufficient, fresh 
young leaves were harvested for DNA isolation. Juvenile 
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growth traits, including seedling height, basal diameter 
and branching number, were measured consecutively for 
3 years at the end of each year. 

Molecular marker analysis 

DNA was extracted from both the megagametophytes 
and young leaves of 96 seedlings by a CTAB method. 
The quality of DNA extracted was tested by measure- 
ment with NanoDrop Sepectrophotomter ND-100 and 
electrophoresis on 1% agarose gels, based on which the 
DNA concentration was calculated. The DNAs of the 96 
samples were used for RAPD, ISSR and AFLP analyses. 

From 100 primers tested in five random samples, a 
total of 16 RAPD primers were found to produce poly- 
morphic markers. Amplifications were carried out in 
20-ul aliquots of the solution containing 12.99 ul H 2 0, 
2ul lOxbuffer, 1.76 ul 25 mM MgCl 2 , 0.3 ul 10 mM solu- 
tion of mixed dNTPs (Sangon Biotech (Shanghai) Co., 
Ltd.), 0.75ul lOuM primer, 2 ul of a 10 ng/ul DNA solu- 
tion, and 0.2 ul 5 U/ul Taq DNA polymerase (TaKaRa 
Biotechnology (Dalian) Co., Ltd.). Amplifications were 
performed in a GeneAmp PCR System 9600 (Applied 
Biosystem, Singapore) programmed for 40 amplification 
cycles (94°C for 0.5 min, 35.2°C for 1.5 min, and 72°C for 
1.5 min) followed by 7-min extension at 72°C. Amplifica- 
tion products were resolved by electrophoresis on a gel 
of 1% agarose run at 110 V for 45 min in lxTBE buffer. 
Band sizes were estimated by comparison to a GeneRuler 
lOObp DNA Ladder (Fermentas). Gels were stained with 
ethidium bromide (EB) and photographed on an Alpha 
Imager (Alpha Innotech Corporation, USA), and the seg- 
regation patterns were scored manually as a band being 
present or absent from a computer printout. The RAPD 
markers were denoted by their primer code followed by 
the numbering of the band (e.g. S30-8). 

An AFLP analysis protocol was optimized based on 
description by Vos et al. (1995). A 20 ul aliquot containing 
250-ng genomic DNA, 2 ul lOxNE-Buffer 2, 0.2 ul 
lOOxBSA, and 12.35ul H 2 0 was digested with the restric- 
tion enzymes viz. 0.15ul of a 20 U/ul £coRI (New England 
Biolabs) and 0.3ul of a 10 U/ul Msel (New England Bio- 
labs) at 37°C for 5 h followed by incubation at 65°C for 
15 min to inactivate the enzyme. The DNA digested frag- 
ments in the 20 ul aliquot were then ligated at 16°C for 
12 hours to 1.4 ul of M-adaptor (50 pM) and 1.0 ul of 
E-adaptor (5 pM) contained in a 5 ul aliquot of O.lul of a 
350 U/ul T4 ligase (TaKaRa Biotechnology (Dalian) Co., 
Ltd.) and 2.5ul 10xT4 ligase buffer followed by a 10-min 
enzyme inactivation at 65°C. The pre-selective amplifica- 
tion cycle was carried out in a GeneAmp PCR System 
9600 (Applied Biosystem, Singapore) programmed for 
25 amplification cycles (94°C for 0.5 min, 56°C for 1 min, 
and 72°C for 1 min) followed by 5-min extension at 72°C 
with one selective nucleotide (E-A, M-C) in an 20 ul 



aliquot containing 2 ul digested fragment solution, 1.5 ul 
M-primer (10 mM), 1.5 ul E-primer (10 mM), 2.0 ul 
lOxPCR buffer, 1.2 ul 25 mM MgCl 2 , 0.3 ul 5 U/ul Taq 
polymerase (TaKaRa Biotechnology (Dalian) Co., Ltd.), 
1.0 ul 10 mM solution of mixed dNTPs (Sangon Biotech 
(Shanghai) Co., Ltd.), and 10.5ul H 2 0. The pre-selective 
amplification product was diluted at 1:30 to be used as a 
template in the following selective amplification. The 
selective amplification was programmed for initially dena- 
turing at 94°C for 3 min., 15 amplification cycles at 94°C 
for 0.5 min, 67°C for 0.5 min, and 72°C for 1 min with a 
lowered annealing temperature of 0.7 each cycle, and 
23 cycles at 94°C for 0.5 min, 56°C for 0.5 min, and 72°C 
for 1 min followed by 7-min extension at 72°C using a pair 
of primers, each of which contained three selective nucleo- 
tides, in a 20 ul aliquot containing 2 ul diluted pre-amplifi- 
cation product, 1.4 ul M-primer (10 mM), 1 ul E-primer 
(10 mM), 2.0 ul lOxPCR buffer, 1.3 ul 25 mM MgCl 2 , 
0.2 ul 5 U/ul Taq polymerase (TaKaRa Biotechnology 
(Dalian) Co., Ltd.), 1.0 ul 10 mM solution of mixed dNTPs 
(Sangon Biotech (Shanghai) Co., Ltd.), and ll.lul H 2 0. 
Totally 41 pairs of primers were used in the experiment. 
The amplification products thus obtained were run on a 
6% polyacrylamide gel for 2.5 h on a sequencer (CBS, 
USA). Gel images were scanned and manually scored for 
the presence or absence of bands. The AFLP bands were 
named after the code of EcoR HMse I primer pairs fol- 
lowed by the number of the band. 

In total, 24 ISSR primers were used in the analysis of 
genomic DNAs. Amplifications were performed in 20-ul 
aliquots of the solution, each of which contained 3 ul of 
a 10 ng/ul template DNA, 2 ul lOxbuffer, 3.25 ul 
25 mM MgCl 2 , 0.5 ul 25 mM solution of each dNTP 
(Sangon Biotech (Shanghai) Co., Ltd.), 0.7ul 10 uM pri- 
mer, 0.2 ul Taq DNA polymerase (TaKaRa Biotechnol- 
ogy (Dalian) Co., Ltd.) and 10.35 ul H 2 0. Amplifications 
were programmed for pre-denaturing at 94°C for 5 min, 
and 35 amplification cycles (94°C for 30 s, denaturing 
for 45 s, and 72°C for 90 s) followed by 7-min extension 
at 72°C, in which the denaturing temperature varied 
among primers. Amplification products were resolved 
by electrophoresis on a gel of 1.6% agarose at 150 V for 
30 min, EB stained and photographed (Alpha Innotech 
Corporation, USA). The segregation patterns of ISSR 
markers were scored and denoted in a way similar to 
those of RAPD makers. 

Map construction and QTL mapping 

Chi-squared tests were performed to check if the markers 
genotyped follow the Mendelian ratios. The markers with 
significant departure from Mendelian segregation ratios 
were excluded from the map construction. A maternal 
linkage map was constructed using MapMaker version 
3.0 [43] following a backcross mapping strategy. The data 
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sets of dominant markers were duplicated to allow the 
detection of repulsion phase between linked markers, 
indicated by R after the names of markers. The "triple 
error detection" and the "error detection" features were 
used to recognize the circumstance in which an event 
was more probably the result of error than of recombina- 
tion. These features avoid map expansion [44]. Linkage 
groups were assigned with thresholds for a minimum 
LOD (logarithm of the odds) score of 6.0 and a maximal 
recombination fraction of 0.30. The markers were 
ordered using the "Order" command iteratively with a 
default LOD of 3.0. The first-sequence order was con- 
firmed using the "Ripple" command permuting five mar- 
kers at a time. All markers that were not ordered in the 
first pass were placed again using the "Try" command. 
Linkage maps were generated with the "map" command 
using the Kosambi mapping function. Maps were drawn 
with the program MapChart 2.1 [45]. 

For megagametophyte-generated markers, we used 
Wu's [18] approach to formulate a mixture likelihood for 
functional mapping of QTLs that control growth traits. 
Functional mapping is a dynamic model that was derived 
to map growth QTLs through integrating growth equa- 
tions into a mapping framework [28-32]. In this study, 
growth traits were measured in the first three years of 
growth. Growth curves that describe juvenile growth are 
defined by an exponential equation [39], i.e., 

g(t) = ae n 

where g(t) is the growth of a trait at time t, a is the 
constant and r is the exponent. Thus, by estimating a 
and r for each QTL genotype, we can test the effect of a 
QTL on growth curves. 

Functional mapping also models the covariance struc- 
ture of growth traits measured at multiple time points 
by a statistical approach. In this study, we implemented 
a first-order autoregressive (AR (1)) model to fit the 
covariance structure by using the variance and the 
between-time correlation. To make the data fitted by 
the stationarity of covariance and correlation assumed 
for the AR (1) model, we log-transformed raw data of 
growth. The log-transformation leads to more parallel 
curves and, thus, increases the feasibility of using the 
analytically advantageous AR (1) model to analyze the 
growth data. 

If different genotypes at a given QTL correspond to 
different trajectories, the QTL must affect the differen- 
tiation of this trait. Therefore, by estimating the curve 
parameters [a, r) that define the trait trajectory of each 
QTL genotype and testing the differences in these para- 
meters among genotypes, we can determine whether a 
dynamic QTL exists and how it affects the formation 
and expression of a trait during development. 
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